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Abstract 

I propose a simple model, based on an analogy to von Neumann artificial viscosity, of turbulent 
diffusion, heat diffusion and viscosity coefficients for use in modeling subgrid turbulent diffusivity 
in multi-phase numerical hydrodynamics and, more generally, in subgrid turbulent viscosity and 
thermal transport. In analogy to the von Neumann artificial viscosity, these coefficients explicitly 
contain the grid size and do not attempt a quantitative model of the unresolved turbulence. In 
order to address the problem that it is often not known a priori when and where a flow will become 
turbulent, the coefficients are set to zero when the flow is not expected to be turbulent on the 
basis of a Richardson/Rayleigh- Taylor stability criterion, in analogy to von Neumann's setting of 
artificial viscosity to zero in expanding flows. 
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I. INTRODUCTION 



Turbulent diffusion in heterogeneous or multi-phase ffows and turbulent viscosity and heat 
transport in all flows are problems of great practical importance in numerical hydrodynamics. 
In many problems the Reynolds number is enormous and if the flow becomes turbulent it is 
not possible to calculate it explicitly on the finer scales of heterogeneity that are still large 
enough that microscopic transport processes (diffusivity, viscosity and conductivity) have 
not smoothed out the variations in the corresponding physical variables. 

When the composition of the fluid is heterogeneous the most important questions usually 
concern the mixing of material of different compositions or in different phases rather 
than the turbulent momentum transfer. A number of sophisticated turbulence models have 
been developed, includ.ng the "K-L- model of Rayleigh-Taylor and Richtmyer-Meshkov 
instability and the more general BHRZ model [3!] (sometimes called the BHR model) that 
implicitly account for all sources driving instability and turbulence, including shear (Kelvin- 
Helmholtz instability), through their derivation from the Euler equations. As powerful as 
these models are, they can be complex to implement because they describe the strength of 
turbulence with parameters that are non-local and history-dependent functions of space and 
time. The success of the very simple von Neumann ^ artificial viscosity, nearly universally 
used in numerical calculations of compressible flow, suggests that a similarly simple model 
of turbulent transport, may be useful. 

II. GENERAL CONSIDERATIONS 

The elementary Kolmogorov turbulence model indicates an eddy turnover time 

tk oc (1) 

for eddy size 1/k. Hence if the flow is turbulent at all small eddies turn over much faster 
than larger ones. Assuming at least one turnover for the largest driving eddies there is likely 
to be rapid mixing down to the inner turbulence scale. 
If the Schmidt number 

Sc = ^, (2) 

where u is the kinematic viscosity and D the diffusivity, is of order unity or less then diffusive 
mixing on scales smaller than the inner turbulence scale will be at least as rapid as the 



turnover of the smallest eddies, and much more rapid than that of the driving eddies. This 
is usually the case because in dilute gases or weakly coupled liquids momentum, thermal 
energy and composition are, to a substantial extent, carried by the same particles (in a 
plasma momentum is almost entirely carried by the ions and energy by the electrons; at 
sufficiently high Reynolds number all of these quantities diffuse on the inner turbulence 
scale much faster than the turnover of the driving eddies). Only if Sc ^ 1 may diffusion 
down to the molecular or atomic scale be slower that turbulent turnover; this may occur when 
the viscosity is large because of strong intermolecular forces, as in viscous liquids (glycerin, 
honey, pitch, the Earth's mantle, etc.) in which momentum diffuses rapidly but chemical 
composition slowly. Such substances rarely undergo turbulent high Reynolds number flows, 
though such flows are, in principle, possible if driven strongly enough. 

In the classical Kolmogorov turbulent cascade the inner scale of turbulence is approx- 
imated by 



where i/ is the kinematic viscosity and e ~ U^/ R is the kinematic (per unit density) energy 
dissipation per unit volume of a flow with characteristic velocity U on an outer scale R. 
This leads to an estimate of the time for microscopic diffusion to produce atomic-scale 
homogeneity 



where the turnover time ttum = R/U. In high Reynolds number flows microscopic (diffusive) 
mixing on length scales from £j„ down to atomic dimensions is a very rapid process unless 
the Sc > 1. 

Turbulent heat transport presents similar problems. If the Prandtl number 



where a is the microscopic heat diffusion coefficient, is of order unity or less then diffusive 
heat transfer on scales smaller than the inner turbulence scale will be at least as rapid as the 
turnover of the smallest eddies, and much more rapid than that of the driving eddies. This 
is usually the case because in dilute gases or weakly coupled liquids momentum and thermal 
energy are, to a substantial extent, carried by the same particles. In liquid metals or ionized 
plasmas the rapid thermal conduction by electrons or rapid radiative transfer lead to Pr <S 1 




(3) 




(4) 




(5) 
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and this inequality holds even more strongly. Just as for turbulent mixing, only in viscous 
liquids in which strong intermolecular forces carry momentum is Pr ^ 1 and the efficacy of 
diffusive heat transport on scales smaller than the inner turbulence scale is problematic; this 
paper does not discuss that regime. 

III. TURBULENCE CRITERION 

The problem is to model the subgrid diffusion in fluids with Sc < 1 when small features 
of the flow cannot feasibly be resolved computationally. This requires a criterion to decide 
when the flow is turbulent as well as a model for turbulent mix. It is often not possible to 
perform an instability analysis because of the complexity and non-stationary character of 
the large scale flow in problems such as inertial confinement fusion. The recognition that a 
complex fiow is turbulent and will possess a turbulent cascade to large wave numbers, much 
less the determination of its characteristics, is non-trivial. This limits the usefulness of even 
the most sophisticated turbulence models. 

Because microscopic diffusion on scales smaller than £j„ is generally much faster than 
the lifetime of the fiow, which is generally 0{tturn), determining whether a fiow is turbulent 
on scales ~ im not explicitly resolved in a numerical calculation is much more important 
than a quantitative estimate of the turbulent diffusion coefficient itself. A common criterion 
for whether a computed or observed fiow is turbulent is "I know it when I see it". This is 
useless in a large numerical simulation in which the structure of the fiow is not known in 
advance. A simple and automatic turbulence criterion is required. Without knowing if a 
fiow is turbulent or not it is impossible even to decide if a turbulence model should be used 
at all, much less to make reasonable estimates of its parameters. 

The model presented here is almost childishly naive compared to the non-oscillatory finite 
volume (NFV) difference schemes js], 0] that have been remarkably successful in modeling 
turbulent flows. So is the von Neumann artiflcial viscosity, yet it remains the basis of most 
numerical calculations of compressible flows. 

The present model uses an on-off switch on turbulent transport based on the Richardson 
stability criterion appropriate to high Reynolds number shear flow (if there is no shear it 
reverts to the Rayleigh- Taylor stability criterion for inviscid miscible fluids). This switch is 
a generalization of one used (j^) in the "K-L" model of turbulence produced by Rayleigh- 
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Taylor and Richtmyer-Meshkov instability. An implicit switch is present in the more general 
BHRZ model jsj of turbulence in which the explicitly calculated velocity field drives or damps 
turbulence according to the Euler equations. 

A switch is not needed in calculations of fiows (such as unstratified pipe fiows) that are 
known a priori to be unstable, and that constitute the usual tests of turbulence models. 
It is essential in more complex heterogeneous fiows that contain both stable and unstable 
regions that are not predictable in advance. The switch is analogous to the switch in the 
von Neumann artificial viscosity that sets it to zero in rarefying regions. 

Sub grid scale (SGS) turbulence models of this type were first proposed by Smagorinsky 
[3] for use in numerical large-eddy simulations (LES) of fiows whose small scale turbulence is 
numerically unresolved. Many SGS models have been developed since his pioneering work 




but none appear to be universally applicable. In this paper I follow von Neumann's 
artificial viscosity model for dissipation in nnresolved shocks, by renonncing any attemf,t 

to model the SGS turbulence in detail or on the basis of fundamental principles, just as [4] 
renounced any attempt to understand the microscopic dissipation mechanisms and structure 
of shocks (which also depend on the detailed physics of the particular fiuids involved, so that 
no single model can be generally valid), and construct a model that instead depends explicitly 
on the numerical properties of the LES, while incorporating a criterion for instability and 
the presence of turbulence and turbulent diffusivity that explicitly depends on the resolved 
properties of the fiow. 



IV. MODEL TURBULENT DIFFUSIVITY 



This problem has some similarity to that of modeling shocks in compressible fiow. In each 
problem the difficulty is posed by subgrid scale dissipative processes that cannot be resolved 
computationally, but are known to be present. The classic solution [4] is to define an "artificial 
viscosity" that is explicitly dependent on the grid size (thus broadening a shock to a width of 
a few resolution elements) and introducing a nonlinear criterion for its presence, in essence an 
on-off switch that turns it on when the fiuid is being compressed and off when being rarefied 
(because in thermodynamically stable systems there are no rarefaction shocks). The resulting 
artificial viscosity is strongly nonlinear and has the unphysical property of a magnitude 
that depends on the scale of the computational grid. It has been widely adopted because 



5 



it accurately reproduces the shock jump conditions, at the price of not reproducing the 
structure of the shock on the scale of the grid resolution or its (much finer and unresolvable) 
actual microscopic structure. It has the same dimensions as real physical viscosity. 

We seek an analogous heuristic artificial diffusivity for multicomponent flows. The tur- 
bulent diffusivity we propose is 

f if > i 

6 p'|e.,|2 ^ 4' 

where a is the local fluid acceleration, A is a characteristic (filter) length, 
is the traceless (anisotropic) symmetrized strain rate tensor |12l . [l3| and 



is its scalar magnitude; summation over repeated indices is implicit. This Dmrb is analogous 
to the Smagorinsky turbulent viscosity model of subscale eddy simulation (see for a 
review) . 



A number of heuristic prescriptions for A have been used [15|, |16|, llTj. Without an 
opportunity to compare to detailed experimental data in complex and difficult to diagnose 
flows such as those encountered in inertial confinement fusion, it is hard to decide which is 
most appropriate for those flows, and the best prescription may vary over space and time. 
We suggest that A be taken to be the smallest dimension of a spatial zone because high 
aspect ratio zones are generally used only when there are large gradients of composition or 
other variables in the direction of the smallest side of a zone that need to be resolved; these 
smallest sides are likely to limit mixing lengths. 

Alternative forms are possible. For example, a tensor Dij could be defined by replacing 
\eij\ by Cij in the expression for Dturh- However, if the turbulent cascade isotropizes the 
turbulence, as generally assumed in the absence of a quantitative more general model, then 
the subgrid Dturb should be a scalar, as in the form above. 

The on-off switch is adapted from the Richardson instability criterion jisl . fiol in strat- 
ified fluids, replacing the stratifying acceleration of gravity g by the component of local 
acceleration a along the gradient of potential density: 

dp 



Vp, (9) 



where the partial derivative is evaluated under thermodynamic conditions appropriate to 
matter displacements. If flows are nearly adiabatic on displacement time scales this is |^ , 
while if they are nearly isothermal (as will be the case if conductive or radiative energy 
transfer to a heat bath is rapid) it is . Stratification can stabilize against either (or 
both) of Rayleigh- Taylor or Kelvin- Helmholtz instabilities. 

This criterion ignores the stabilizing effect of viscosity, which is small at high Reynolds 
numbers and not easy to calculate when both shear and stratification affect stability. At 
lower Reynolds numbers (defined by zone size), viscosity is significant but calculated directly 
by a hydrodynamics code. When this is the case subgrid scale turbulence does not occur 
and, provided the Schmidt number is 0{1) molecular diffusivity is rapid and no turbulent 
diffusivity model is required. 

The factor of 1/6 is introduced to allow for three spatial dimensions (in analogy with the 
factor 1/3 in the elementary kinetic theory result D = ui/3, where i is the mean free path), 
and for the double counting resulting from symmetrization of Cij. Because the proposed 
form is phenomenological, if fitting data are available it may be appropriate to introduce 
an additional dimensionless multiplicative factor chosen to match those data. Such "data" 
may be obtained from a fine scale but idealized 3-D numerical hydrodynamic calculation of 
turbulent mix, for example in an unstratified isochoric (constant and uniform density) fluid. 

When this model is applied to an Eulerian calculation it is tacitly assumed that an 
interface-preserving algorithm is used to prevent rapid numerical diffusion of composition 
that may be entirely spurious (for example, if the stratification is stable). 

V. DISCUSSION 

Turbulent mixing is a problem of long standing, and many models have been introduced. 
Typically they involve several free parameters, and these may be chosen to give excellent 
fits in regions of parameter space for which appropriate data are available. For example, the 
widely used BHR modeljsl, [2O] contains ten fitting parameters. These have been chosen to fit 
a classic shock tube experiment jsil . With so many parameters available the fit is excellent, 
but that does not establish its validity when extrapolated outside the range of conditions to 
which they were fitted. For many important problems no data exist in the parameter range 
of interest. A simple model may be more robust than a sophisticated one. 
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Any model, and especially one as oversimplified as that proposed here, needs to be tested. 
Quantitative experimental tests are few because of the difficulty of measuring (and even 
defining in an experimentally accessible manner) effective diffusion coefficients for complex 
heterogeneous flows. Numerical simulation may be a preferable approach, because there is 
no difficulty in defining spatial distributions of composition, and even of fitting an effective 
D to the evolution of their means. An unstable velocity field is set up on a fine mesh, and 
the mixing calculated by the code is compared to that produced by the model applied to a 
much coarser description of the mean velocity field. Both the validity of the model and the 
optimal value of any multiplicative coefficient may be determined from this comparison. 

Subgrid turbulent viscosity and heat transport present problems analogous to that of 
subgrid turbulent mixing. Sophisticated models exist, but also contain several parameters 
(essentially closure parameters for the moment equations). It may be a reasonable rough 
approximation to use the model proposed here for Dturb and to take i/turb — ^CturbDturb, 
estimating Scturb from numerical simulations like those used to test the model for Dturb, 
or simply taking Uturb = Dturb in the absence of such information. Similarly, it may be 
reasonable to take aturb = Dturb- 

I thank G. Dimonte, P. E. Dimotakis, R. A. Gore, L. G. Margolin, D. I. Meiron and 
G. B. Zimmerman for useful discussions. 
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